Impact of the 2008 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_W$$\end{document}MW 7.9 Great Wenchuan earthquake on South China microplate motion

Tectonic plate motions drive the earthquake cycle, as they result in the slow accrual and sudden release of energy along plate boundaries. Steadiness of plate motions over the earthquake cycle is a central tenet of the plate tectonics theory and has long been a main pillar in models of earthquake genesis, or of plate-margins seismic potential inferred from slip-deficit estimates. The advent of geodesy in the geosciences and the availability of multi-year-long series of position measurements permit tracking the motions of tectonic plates from before to after the time of significant seismic events that occur along their margins. Here, we present evidence that large earthquakes are capable of modifying the motions of entire microplates. We use high precision Global Navigation Satellite System (GNSS) position time-series covering the periods 2001–2004 and 2014–2017 to demonstrate that, contrary to the tenet above, the South China microplate motion changed after the 2008 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_W$$\end{document}MW 7.9 Great Wenchuan earthquake. The GNSS data and associated uncertainties indicate a plate motion slowdown of up to 20% that is beyond the possible impact of data noise and is thus tectonically meaningful. We use quantitative models of torque balance to show that generating this kinematic change requires a force upon the South China microplate compatible with that imparted by the Great Wenchuan earthquake of 2008. The existence of a kinematic signal linked to the earthquake cycle that impacts an entire microplate might offer an additional, novel perspective to assessing the hazards of earthquake-prone tectonic regions.

The earthquake cycle is a process whereby relative motions between adjacent tectonic plates determine the slow, decade to century-long buildup of stress along tectonic plate margins followed by sudden stress release through earthquakes 1 .The great majority of large, destructive earthquakes occurring within Earth's crust can be studied in this context.A central assumption in the geosciences is that the plate motions that fuel the earthquake cycle are insensitive to the temporal variations of crustal fault stresses associated with such a cycle.This assumption is the foundation of most numerical models of spontaneous earthquake generation that implement rate-and statedependent frictional laws [2][3][4][5] , or models of the seismic potential of global plate margins that employ estimates of slip-deficit [6][7][8][9][10][11] .Geodetic observations, in particular the Global Navigation Satellite System (GNSS), provide the opportunity to measure contemporary plate motions 12,13 .The use of GNSS position time-series spanning periods shorter than the typical duration of earthquake cycles 1 allows resolving plate motions at different moments during such cycles.Recent analyses of GNSS data in the Mediterranean and Andean settings [14][15][16][17] indicate non-steadiness of plate motions during the earthquake cycle and thus call for testing the tenet above in other tectonic settings.
The South China (SC) microplate is an independently-moving continental lithospheric block located adjacent to the Eurasia plate, between the Tibetan Plateau and the trench where the Philippine Sea plate subducts beneath Eurasia 18 (Fig. 1).The present-day SC motion is directed towards the East/Southeast relative to the Eurasian plate, as evidenced from contemporary GNSS measurements 19,20 , and is driven chiefly by the extrusion of Tibet 21 .Analyses of GNSS data 22 indicate that most of SC is characterised by little internal deformation, with baseline linear strain-rate estimates around 1.5×10 −9 yr −1 .In contrast, the continental crust surrounding SC features strain-rates between 3 and 10 times larger 13 , which makes a precise identification of portions of the SC tectonic margins not straightforward 18,23 .On 12 May 2008, part of the western margin of SC ruptured with a M W 7.9 seismic event that remains one of the strongest recorded in China since the middle of the 19 th century and is referred to as the Great Wenchuan earthquake.The earthquake occurred on the Longmen Shan Fault system and .Velocities observed at continuously-recording (solid arrows) and campaign (empty arrows) GNSS stations used in this study to infer Euler vectors for the South China (SC, in black) microplate motion relative to Eurasia during the periods from January 2002 to December 2004 (in red), and from January 2015 to December 2017 (in blue).All GNSS stations are labeled: italics is used for campaign GNSS stations available during the earlier time period, bold for continuous GNSS stations, and bold-italics for the four continuously recording sites available for both the earlier and later time periods.The black beach ball shows the epicentre and USGS focalmechanism solution of the 2008 M W 7.9 Great Wenchuan earthquake.Turquoise dots are aftershocks.This figure was created utilising the Generic Mapping Tools version 6 (https://www.generic-mapping-tools.org/) 64 .
analyse GNSS position time-series collected at sites of the Crustal Movement Observation Network Of China (CMONOC) 22,34 over two 3-and 3.5-year-long periods, before and after the 2008 Wenchuan earthquake.We calculate GNSS velocities that we use to constrain statistically significant temporal changes of the Euler vector for the SC microplate motion from before to after the 2008 earthquake.Lastly, we estimate the force required to generate the inferred kinematic change.We find that the estimated force is comparable with that imparted to SC by the 2008 Great Wenchuan earthquake coseismic stress drop.

GNSS data
The motion of SC can be expressed using Euler vectors 35 , which is an approach that assumes minimal internal plate deformation-smaller than the uncertainty of velocity measurements constraining the Euler vector.We constrain SC motion from GNSS sites belonging to the CMONOC 22,34 that are located within the SC margins, in most cases at least 100 km from the plate boundary.The CMONOC was established in 1998 and significantly developed after 2010.Specifically, we focus on data recorded at the available GNSS sites over two time periods, one before and one after the 2008 Great Wenchuan earthquake.The time period before the event is from 1 January 2002 to 25 December 2004 (in the following referred to as earlier period), and the time period after the event is from 1 January 2015 to 25 December 2017 (in the following referred to as later period).The 3-year-long periods are selected using the following criteria: (i) they are longer than 2.5 years, which is generally considered the minimum standard to obtain tectonically-representative motions 36 .(ii) Following previous studies 14 , the later time period begins at least one Maxwell time-interval of the sub-SC asthenosphere (see "Methods" and Supplementary Information) after the earthquake.The occurrence of the 2008 coseismic stress drop means that a force acting upon SC (one that has grown during the interseismic phase over the area that later ruptured) is suddenly removed.This causes a readjustment of the motion of the SC lithosphere that, in turn, exerts a sudden stress at the lithosphere-asthenosphere interface.The asthenosphere responds to that in a viscoelastic fashion 1,37 -meaning that there is an immediate (after the 2008 earthquake), predominantly-elastic response that transitions through time into a predominantly-viscous response.The time it takes for this transition to occur is of the order of the Maxwell time interval.Previous studies 38,39 have shown that the postseismic deformation in SC following the 2008 event is negligible.Nonetheless, we take a cautious approach by choosing a later time period that begins after one Maxwell time interval from the 2008 earthquake.This is to ensures that GNSS sites measure velocities at a time when most of the already negligible post-seismic elastic response of the viscoelastic asthenosphere beneath SC has elapsed, and asthenospheric flow is characterised by a purely-viscous component.In other words, we minimise the possibility that the GNSS sites record any of the post-seismic elastic part of the response 37 .(iii) By the same logic, the selection of earlier/later periods aims at minimising any possible postseismic effects of other large earthquakes occurred around SC, such as the 1999 Chi-Chi, the 2004 Sumatra, and the 2011 Tohoku earthquakes.Several independent studies indicate that the postseismic-relaxation signals associated with these events over the earlier/later periods become negligible if one moves as far away from the epicenters as SC is [39][40][41] .The fact that de-trended position time-series (available in the Supplementary Information) do not exhibit monotonic trends confirms those inferences.(iv) Lastly, the basis of our working hypothesis is that microplate motions might be impacted by temporal stress variations associated with any seismic cycle that occurs along their margins and that precedes large earthquakes.This hypothesis implies the possibility that the SC motion recorded at GNSS sites is impacted by other stress variations along the SC margins that are not part of the 2008 Great Wenchuan earthquake cycle, but are associated with the interseismic phase preceding large earthquakes yet to occur.A previous study 42 used synthetic models of plate dynamics to determine that while in principle the seismic cycle associated with any earthquake occurring along the margins of a microplate has an impact onto its motion, only the impact of large earthquakes is indeed sufficiently large to change the microplate motion by more than the typical uncertainty associated with GNSS inferences.This suggests that contemporary GNSS motions of plates should not be expected to remain consistent to their long-term (i.e., derived from paleomagnetic observations) counterpart, but rather experience short-term changes in response to seismic cycles going on along their margins.For the reasons above, we select the later period to be relatively close in time to the 2008 coseismic stress drop, even if position time-series recorded at sites of CMONOC extend up to today.Our choices for the GNSS observation time periods assumes that interseismic stress might grow substantially in the latest interseismic phase 43 .Thus, with our GNSS observation time period choices, we aim to minimise potential superimposed effects of ongoing interseismic stress variations caused by unknown seismic cycles that have not yet reached the coseismic phase.
During the earlier time period-from 1 January 2002 to 25 December 2004-4 continuously-recording GNSS stations belonging to CMONOC were deployed in SC.During the later period-from 1 January 2015 to 31 December 2017-there were 29 continuously-recording GNSS stations in SC-including the 4 continuous stations available for the earlier period.These stations represent the entirety of continuous GNSS data ever collected within SC over these periods.Given the marked difference in number of continuous GNSS stations available for the earlier and later periods, we elect to complement the earlier period with CMONOC campaign GNSS data.We consider CMONOC campaign sites that have been recording between 1999 and 2004, fit a linear model, and are located at least 100 km away from the SC tectonic margins to ensure they are not affected by plate boundary deformation 13 .For the earlier period, CMONOC provides 17 campaign GNSS sites that meet such criteria.We process position time-series from these sites over the earlier and later time periods (see "Methods") to obtain velocities and associated uncertainties (Fig. 1).All position time-series and velocities are linked in the Additional information section.
The continuous GNSS stations LUZH, GUAN, XIAM, and WUHN are common to both the earlier and later periods (see sites with full red and blue arrows in Fig. 1).In addition, 4  www.nature.com/scientificreports/JB25) that recorded during the earlier period are located less than 30 km away from some of the continuous stations (GDSG, ZJWZ, SCNC, and GZGY) that recorded during the later period.Velocities at these sites provide a basis for a preliminary assessment of our hypothesis: we compare GNSS velocities observed during the earlier and later periods at the continuous stations LUZH, GUAN, XIAM, and WUHN, as well as at the 4 pairs of closely-located campaign/continuous sites (F037-GDSG, JB19-ZJWZ, JB24-SCNC, and JB25-GZGY)-see Supplementary Fig. 1 for the locations of these sites.We do so by generating ensembles of 10 6 samples of East/North components of velocity using the nominal values and associated standard deviations.From these ensembles, we generate distributions of the temporal change of velocity components calculated by subtracting the ensemble of velocities during the earlier period from that of the later period (Fig. 2).Had there been no temporal change of velocity, the distributions would be centred at the origin-i.e., (0,0) in the 2D space of East/North velocity components.However, 5 out of 8 comparisons exhibit a slowdown of motion that is consistent with our hypothesis (Fig. 2).We note that the remaining 3 distributions, which exhibit microplate motion speedup, are systematically located in the part of SC featuring the longest Maxwell time-interval.This observation suggests that the velocities recorded at these sites during the later period (i.e., the one after the 2008 Wenchuan earthquake) may continue to be impacted by residual postseismic deformation.We note that the velocity comparison at site XIAM (Fig. 2d) does not support significant temporal change of motion at stricter levels of confidence, but we follow Gordon 35 and test our hypothesis based on the analysis of SC Euler vectors, rather than any one single-site observation.
On the basis of the arguments above, we proceed as follows: we initially infer Euler vectors for SC microplate motion over the 3-year-long earlier and later periods starting from velocities from all available continuous GNSS stations.We compare Euler vectors and associated uncertainties, which propagate from standard deviations of site velocities.Next, we evaluate the robustness of any inferred temporal change using two independent statistical tests that weigh the number of observational constraints against the degrees of freedom of the inversion problem: the F-ratio and noise tests.Furthermore, we re-calculate Euler vectors for SC microplate motion also using velocities observed at continuous sites over 3.5-year long periods.Lastly, we repeat the Euler-vector inference in two additional cases: using position time-series from both continuously-recording and campaign GNSS sites over 3-year-long periods, as well as excluding from the analyses those sites located in the northwestern portion of SC featuring a higher Maxwell time-interval.Altogether, these analyses provide a basis for accepting or rejecting any evidence of temporal change with more confidence than that based on sole propagation of site-velocity uncertainties.

Euler vectors for the South China microplate motion
We constrain Euler vectors and associated covariances that represent the motion of SC relative to Eurasia before and after the Great Wenchuan earthquake as follows: for each time period, we start by generating 10 6 samples of For the latter ones, data from the first site (campaign) in the label cover the earlier period, while data from the second site (continuously-recording) cover the later period.Distributions are obtained from the difference between ensembles of site velocities sampled from the nominal values and standard deviations of the East and North velocity components.Red dots show the origin (0,0), which would correspond to the centre of the distribution in case of no temporal change of site velocity.Contours show 20%, 68%, and 95% most-recurrent portions of the distributions.See Supplementary Information for maps of sites and site-pairs.site velocities from the nominal solutions and their uncertainties.We then use these samples to solve the wellknown inverse problem based on minimisation of the sum of squared misfits 44 (see "Methods"), and obtain an ensemble of Euler vectors describing the SC motion.From these ensembles, we calculate mean values and covariances of the Cartesian components across the ensembles.After inverting for the best-fitting Euler vector, we check if all site velocity-residuals (i.e., the vectorial difference between observed site velocities and velocities predicted by the Euler vector at each site) are smaller than their 1-sigma uncertainties.If not, we iteratively repeat the inversion, each time removing, from the set of utilised sites, the one site yielding the largest unacceptable (i.e., larger than 2 standard deviations) residual in the previous iteration.The process stops when residuals of all utilised sites are acceptable.This approach is deliberate, and predicated on the following: realistically, the velocity sampled at each GNSS site is the SC motion, plus some (ideally negligible) contribution of other processes unrelated to the SC motion-whether geological in origin or not.The inversion scheme for Euler vectors weighs all utilised sites equally.On the one hand, the larger the number of utilised sites is, the more precise (i.e., smaller covariance) the Euler-vector inference is.On the other hand, if a specific site features a significant velocity-component unrelated to the SC motion, then such a component (i) impacts negatively onto the accuracy of the inferred Euler vector, and (ii) makes the velocity-residual of the specific site larger than the uncertainty.Our approach is to avoid the latter type of sites by iteratively removing them, and to to verify a-posteriori via noise tests (see below) that any undetected presence of non-plate-motion velocity-components in the pool of remaining sites is nonetheless innocuous to the working hypothesis.This approach carries a gain in terms of accuracy.The associated cost in terms of precision is in our opinion acceptable in light of the noise tests performed a-posteriori.In all cases related to the earlier period, residuals are acceptable when using all sites in the starting set of GNSS velocities, thus there is no need to remove any site in the Euler-vector inversion.In cases related to the later period, we iteratively remove up to 3 sites from the initial set of velocities (Table 1, Fig. 3 and Supplementary Information).The fact that no, or very few, sites need to be excluded in order to obtain acceptable residuals upon inversion for the best-fitting Euler vectors confirms a posteriori that the utilised networks of GNSS sites, indeed, record mostly the SC microplate motion.
Fig. 4A shows a comparison of the ensembles of Euler vectors obtained for the earlier and later 3-year-long periods when using only continuous GNSS stations.Similar figures illustrating results for the other cases explored are reported in the Supplementary Information.The distributions of Euler poles and angular velocities for the two time periods indicate a clear difference between the Euler vectors of the earlier period and the later period.The mean angular velocity for the later time period (in blue) falls outside the 68% confidence interval around the mean angular velocity of the earlier period (in red).Differences between distributions of earlier and later Euler vectors remain evident-although they do not necessarily follow the same patterns shown in Fig. 4A-also in the other cases explored (see Supplementary Information).

Tests for the robustness of the inferred Euler-vector temporal change
Given two sets of GNSS velocities, the widely-employed F-ratio test 45 assesses the probability that one single Euler vector would fit both sets to the same level or better than two distinct Euler vectors.We perform the F-ratio test on the pair of Euler vectors in Fig. 4A, and find there is less than 1.5% probability that a single Euler vector (i.e., one that has not changed through time) would fit the two sets of GNSS site-velocities as well as two distinct Euler vectors.To further test our hypothesis, we perform an additional set of statistical tests where we (i) assume-counter to our hypothesis-that there has been no temporal change of SC Euler vector, and (ii) assess the probability that the Euler-vector change evidenced by GNSS data is only apparent and owes to random, sitecorrelated data noise.Specifically, we perform semi-synthetic tests 14 where we utilise computer-generated noise values and apply them to the SC motion inferred from actual GNSS velocities.In these tests, noise is intended as any site-velocity component that is unrelated to the SC microplate motion.Noise arises from other processes that may or may not be geological in origin.
We start by assuming that the SC Euler vector has not changed from before to after the 2008 earthquake, but has remained steady and equal to the Euler vector inferred from GNSS velocities for the period before the 2008 earthquake (labeled as B3.0P in Table 1).This assumption means that we assign to B3.0P the role of the true Euler vector and to the GNSS sites utilised for the period after the 2008 earthquake the role of noisy sites.By implication, the discrepancy between B3.0P and the Euler vector obtained for the post-Wenchuan period-A3.0 in this case-is assumed to be solely the result of data noise.Next, we use Euler vector B3.0P to calculate surface velocities at each of the GNSS sites used for the equally-long (i.e., 3 years) period after the 2008 earthquake.These are the velocities that GNSS sites would have recorded during the later period under the assumptions above and in the absence of any noise.We then add to each of these site velocities an ensemble of 10 6 values of noise that have been randomly sampled from Gaussian distributions whose standard deviation is a free parameter in range from 10 −2 to 8 mm/yr (referred to as noise standard deviation).Following previous studies 46 , we assume noise is correlated across sites and test correlation coefficient values of 0.1 and 0.9.For each value of noise standard deviation, we generate ensembles of 10 6 noisy site velocities and calculate from these ensembles 10 6 realisations of noisy Euler vectors.Next, we quantitatively assess how close the 10 6 samples of noisy Euler vectors are to A3.0: we count what fraction of the ensemble of noisy Euler vectors (i.e., how many of the 10 6 samples) falls inside a contour/ interval (referred to as target range) where the most-recurrent 68% percent of the Euler pole/angular velocity of the A3.0 ensemble are located.Considering the size of the ensembles, such a fraction arguably represents a reasonable estimate of the probability that noise in the data utilised for the period after the 2008 earthquake is responsible for making the inferred Euler-vector appear different from that obtained using GNSS velocities for the time period before the earthquake.We repeat the same test, this time assigning the role of true Euler vector to A3.0, and the role of noisy sites to those utilised for the period before the 2008 earthquake.Results show that such probability remains low (Fig. 4B, statistics associated with these tests in the other cases explored here are reported in the Supplementary Information).As part of these assessments (i.e., under the assumption that no actual, tectonically-meaningful change in SC microplate motion has occurred), we also perform F-ratio tests to sample the probability that, despite the presence of data noise, velocities observed before and after the 2008 earthquake can be fit by one single Euler vector (i.e., no temporal change of Euler vector) as well as by a set of two distinct Euler vectors.This serves to investigate the level of F-ratio-test probability that one can expect if there was no real change in the SC motion, and the only factor causing the inference of an apparent Euler-vector change was noise.The obtained probabilities exceed 85% (inset in Fig. 4B and Supplementary Information).These results suggest that it is unlikely that there has been no change in the SC motion following the 2008 Great Wenchuan event.The probability that (i) a single Euler vector is adequate to fit GNSS velocities recorded before and after the 2008 earthquake (i.e., that the whole-SC motion has not changed from before to after 2008), and (ii) the apparent change inferred when two distinct Euler vectors are fit to data (Fig. 4A) owes   1) and after (A3.0 in Table 1) the 2008 Great Wenchuan earthquake, inferred from continuously-recording GNSS velocities (Fig. 1).Map shows the positions of the average Euler poles inferred for the two time periods, together with the contours where the most-recurrent 68% (dashed lines) and 95% (solid lines) of the sampled Euler poles fall.Continents are in light grey.SC margins are in dark grey.Inset shows the distributions of sampled angular velocities associated with the Euler poles.(b) Estimate of the probability that random, cross-site correlated noise in GNSS velocities is the only factor determining the difference between inferred Euler vectors for the earlier and later periods.Inset shows the histogram of the obtained F-ratio-test probabilities in noise tests (see main text for details).Part of this figure was created utilising the Generic Mapping Tools version 6 (https://www.genericmapping-tools.org/) 64 www.nature.com/scientificreports/exclusively to data noise is the product of the probabilities obtained in the F-ratio and noise tests, which is at most 0.02% (Fig. 4B).Taken together, these statistical analyses support the notion of a real, tectonically-meaningful temporal change of the SC microplate Euler vector.Using the Euler vectors for the earlier and later periods, we calculate surface velocities at locations that are arbitrarily selected and uniformly spaced (Fig. 5).This illustrates how the Euler-vector temporal change inferred from GNSS data corresponds to a ∼20% slowdown of the East/Southeast-directed SC microplate motion.Sur- face velocities calculated from Euler vectors in the other cases explored here are reported in the Supplementary Information, and systematically illustrate a slowdown larger than that permitted by Euler-vector uncertainties.

Torque-variations upon the South China microplate
Any temporal change of the SC Euler vector must arise from a change in the torques acting upon SC.Parameterising the torque-balance equation for tectonic plates at two moments in time, and then taking their difference, yields a relationship that links changes in torque acting upon a plate to the associated change of the Euler vector via a linear operator (see "Methods").The linear operator accounts for the geometry of the plate as well as the viscous resistance of the asthenosphere at its base 14,42,47 .We use the MYRIAM software 48 to apply such parameterisation to the temporal change of the SC Euler vector in order to constrain the torque-variation needed upon SC to generate a 20% slowdown of its motion.Specifically, we utilise the ensembles of 10 6 samples of the Euler vectors describing the motions of SC before and after the 2008 earthquake to generate ensembles of torque-variation vectors.Furthermore, we repeat such calculation in a second scenario, where we assume that the plate-tectonic approximation (i.e., the use of Euler vectors to describe the motion anywhere within a plate) progressively fails within 100 km of the SC margins (see "Methods").This assumption means that crustal deformation is deemed non-negligible within such a 100-km-wide buffer zone around the SC margins.The distribution of the ensemble of torque-variation vectors can be visualised by focusing on the magnitude of the torque-variation and the geographical location where the torque-variation axis intersects Earth's surface-also referred to as the torque-variation pole (Fig. 6).We sample the distribution of torque-variation required upon SC assuming radial and lateral heterogeneities in the viscosity of the asthenosphere.These heterogeneities are mapped from global surface-wave seismic tomography 49 and scaled assuming end-member values for the global average viscosity of the asthenosphere equal to 1 × 10 19 and 3 × 10 19 Pa s 50 (see "Methods").
We compare the distribution of torque-variation required to act upon SC to generate the observed motion slowdown with the torque-variation imparted to SC by the coseismic stress drop of the Great Wenchuan Surface velocities predicted from Euler vector B3.0P for the SC/EU motion during the period from 1 January 2002 to 25 December 2004.
Surface velocities predicted from Euler vector A3.0 for the SC/EU motion during the period from 1 January 2015 to 25 December 2017. .
earthquake.The latter torque-variation is calculated from models of either the finite-fault slip, which are constrained by teleseismic waves observations [25][26][27][28] , or of the focal mechanism [29][30][31][32] .In calculations relying on finitefault solutions, we assume a linear relationship between coseismic slip and stress drop over the rupture area 51 .
In calculations relying on focal-mechanism solutions, we assume a uniform stress drop over the rupture area (see "Methods").These different approaches allow for testing the working hypothesis despite the large uncertainty on the actual value of the coseismic stress-drop and its lateral distribution over the rupture area.We test maximum stress-drop values in range from 0 to 2 and from 4 to 6 MPa, which are in line with global statistics 52 .
Both the pole and magnitude of torque-variation required to act upon SC are compatible with what is imparted to SC by the rupture of the Great Wenchuan earthquake.This compatibility indicates that the force required to generate the observed slowdown of SC motion and the one provided by the 2008 earthquake are oriented in the same way.The compatibility also indicates that the two forces are of similar magnitude, regardless of the tradeoffs between the controlling parameters that inevitably impact the torque calculations.These inferences also hold when using a hybrid set of campaign and permanent sites for calculating the Euler vector for the earlier period (Fig. 7); when excluding from the Euler-vector inference GNSS sites located in the SC region featuring a high Maxwell time-interval (Fig. 8); when using Euler vectors inferred from GNSS velocities calculated over a 3.5-year-long period (Fig. 9); as well as when accounting for possible non-negligible deformation near the SC margins (dotted contours/lines in Figs. 6, 7, 8, 9).The existence of a whole-plate kinematic signal associated with Figure 6.Comparison of the torque variation required upon SC to generate the plate-motion change evidenced by GNSS data (Euler vectors labeled as B3.0P and A3.0 in Table 1), and the torque variation provided to SC by the coseismic stress release of the 2008 Great Wenchuan earthquake.Acronyms FFS/FMS preceding the source reference stand for finite-fault/focal-mechanism solution, respectively.Map shows the regions where poles of the sampled torque-variation fall.Inset shows the distributions of the magnitude of the sampled torque-variations required upon SC, and compares them to ranges of torque-variation magnitudes provided by coseismic stress release.In hatched is the SC microplate boundary.
the occurrence of large earthquakes, combined with the ever-growing body of data provided by the network of continuously-recording GNSS sites 53 might provide an additional perspective to efforts aimed at assessing the seismic potential of tectonic settings.

Methods
Some of the methods presented below are reprised from previous studies (see references).They are reported here for completeness and reproducibility.

Viscosity and Maxwell time interval of the asthenosphere
We resort to previous models of long-wavelength glacial rebound data 54,55 that constrain the relationship between viscosity and thickness of the asthenosphere.Specifically, they found that the cube of the asthenosphere thickness ( H a ) is proportional to the ratio of asthenosphere viscosity ( µ a ) over viscosity of the lower part of the upper mantle ( µ M ).That is, H a = a • (µ a /µ M ) 1/3 , with a = 4.73 • 10 5 m .Setting the global average of µ a to 5 • 10 19 Pa s and that of µ M to 1.  (15), and ( 17) in their study).We also utilize PM_v2_2012 to map temperature-dependent lateral variations of the asthenosphere Young's modulus, whose global average within the asthenosphere is set to 1.5 × 10 11 Pa 56 .Asthenosphere viscosity and Young's modulus are used to infer lateral variations of the char- acteristic Maxwell time-interval of the asthenosphere.See Supplementary Information for figures illustrating the vertically-averaged asthenosphere viscosity, Young's modulus, and Maxwell time-interval underneath the South China microplate.

GNSS data processing
The GNSS data we processed are from the Crustal Movement Observation Network of China (CMONOC), operated by the China Earthquake Administration.We processed the CMONOC data together with data from the surrounding International GNSS Service (IGS) network to obtain loosely constrained daily solutions using the GAMIT/GLOBK software version 10.7 57 (RELAX configuration), each with associated covariance matrices.We employed the latest antenna phase center models and IGS final satellite orbits in the International Terrestrial Reference Frame 2014 (ITRF2014) reference frame 58 .We utilized the empirical GPT2 59 to determine a prior zenith tropospheric delay and Vienna Mapping Function 1 (VMF1) models for mapping the slant delays along each receiver-satellite path to the zenith.The local loosely constrained daily solutions and their covariances were then combined with global IGS solutions from Scripps Orbit and Permanent Array Center (SOPAC), and then fixed the combined daily solution into the 2014 International Terrestrial Reference Frame 58 .Initially, the combined time-series are edited: time-series points with sigma greater than 2 cm in either horizontal component or 4 cm in the vertical component are excluded preliminarily.Time-series points are deemed outliers and thus removed iteratively if the ratio of Res/(Sigma*NRMS) is greater than 2.5, where Res is the residual time-series point (after removing from the raw time-series a linear rate and offsets at breaks plus any user defined parameters), Sigma is the observed uncertainty of the time-series point, NRMS is the normalised root-mean-square of the residual time-series.The cleaned time-series are then used to derive site velocities in a Kalman filter estimator of trajectory models.Specifically, the trajectory models we employ include known discontinuities (resulting from instrument changes or other sources), a linear term (whose slope is the site velocity), as well as an annual term.Plots of the detrended time-series (i.e., observed series minus trajectory model) are reported in the Supplementary Repository associated with this study (see Additional Information section).For the later period, they show negligible postseismic deformation, in line with previous studies 38,39 .Realistic uncertainties of the velocities of the sites Figure 8. Same as Fig. 6, but for Euler vectors B3.0PC* and A3.0* in Table 1.
Figure 9. Same as Fig. 6, but for Euler vectors B3.5P and A3.5 in Table 1.where the right-hand side of equation ( 12) is known from observations.As observations v i typically come with uncertainties, one can first use these to generate an ensemble of N ω samples of each v i , and then apply the inver- sion scheme above to generate an equally large ensemble of ω .Realistic covariances associated with ω can also be calculated from such ensemble, provided that N ω is sufficiently large.In this study, we use N ω = 10 6 and assume no correlation between the East-and North-directed components of v i .A set of Matlab functions that implement the inversion illustrated above is released as part of this study (see Supplementary Information).Table 1 reports Euler vectors and associated covariances inferred via the method above.

Torque variation required to generate plate-motion changes
Previous studies 42,47 obtained an analytical equation linking the torque variation M experienced by a tectonic plate of basal area S to the resulting temporal change of Euler vector.Writing two equations that express the torque balance of a tectonic plate at two distinct points in time, t 1 and t 2 , and then taking the difference between them, so that only the terms that have actually changed through time remain in the resulting equation, results in where µ a and H a are viscosity and thickness of the asthenosphere, v p (r, t) is the plate motion at position r and time t 42,47 , while dS is the infinitesimal element of plate basal surface S centred around position r .The plate- tectonic description means that v p (r, t) = ω(t) × r , which would make Eq.( 13) become where �ω = ω(t 2 ) − ω(t 1 ) .The possibility that deformation near plate margins may be non-negligible would imply that surface velocities and their temporal changes near plate margins would not strictly follow the Eulervector description.In order to account for this, one may introduce the option of locally dampening temporal kinematic changes implied by Euler-vector changes.To very first order, this notion can be implemented by means of a location-dependent dampening factor D r ≤ 1 .D r = 1 if position r is sufficiently away from plate margins, while D r decreases from 1 to 0 as r gets closer to plate margins.Thus, Eq. ( 14) becomes The fact that the Euler-vector change �ω is independent of r means that the integral may be written as the result of a linear map where the operator P links the vectorial space of Euler-vector changes that a tectonic plate could experience to that of torque variations possibly acting upon it.Therefore, Equation ( 17) holds when the time periods over which M and �ω are accomplished are (i) short enough (in a geological sense) that the plate has not changed shape, but also (ii) long enough that the elastic component of sublithospheric motion within the viscoelastic asthenosphere-which affects velocities of the overlying plate 37 -has decayed to a significant extent.These conditions indeed apply to SC for Euler vectors calculated using GNSS data collected between 2001 and 2004, and between 2014 and 2017, since these two periods are (9) �M = P�ω more than 2-3 Maxwell time-intervals apart from each other (see Supplementary Information).We draw 10 6 samples of the SC Euler vectors (Table 1) for periods before and after the Great Wenchuan earthquake, and use them to build ensembles of geodetically-constrained SC Euler-vector changes that we then map, through equation (17), into ensembles of torque variations needed to generate the SC temporal change of motion.In doing so, we implement the findings of previous studies of long-wavelength glacial rebound data 54,55 , which constrain the cube of the asthenosphere thickness to be proportional to the viscosity contrast between the asthenosphere and the upper mantle.In addition, we also implement asthenosphere lateral variations of viscosity in line with the lateral temperature variations mapped by the tomography model PM_v2_2012 49 .Furthermore, we explore the impact of two distinct scenarios for the dampening factor D r : in one scenario, we set D r = 1 everywhere inside S.This means assuming SC as perfectly rigid.In the second scenario, we assume a 100-km-wide internal buffer zone along the SC margins where deformation is non-negligible.Inside this zone, D r decreases linearly from 1 to 0 as r gets closer to the SC margins, and becomes equal to 0 on the SC margins.The entries of the SC-specific operators P when the viscosity of the upper mantle is set to 1.5 • 10 21 Pa s and µ a is set variably to 1 • 10 19 and3 × 10 19 Pa s 50,63 are reported in the Supplementary Information.We visualise the distribution of M by plotting the histogram of its magnitude | M| and the contour containing 95% of the geographical positions where the sampled torque-variation directions �M/|�M| intersect Earth's surface-also referred to as torque- variation poles.All the above is performed via the MYRIAM open-source software 48 .

Torque variation from co-seismic stress drop
Following previous studies, we calculate the torque variation imparted to SC by the 2008 Great Wenchuan earthquake rupture ( M e ) using two methods, one utilising finite-fault slip solutions, the other utilising focal- mechanism solutions.In the former case, we follow the method of Martin de Blas et al. (2022) 14 and proceed by integrating the coseismic stress drop occurred during the earthquake over the finite rupture area.We assume that the stress drop over the rupture area scales linearly with the amount of coseismic slip 51 , which has been estimated through teleseismic wave observations [25][26][27][28] and is available in the form of finite-fault rupture solutions at http:// equake-rc.info/ srcmod/ 27 .In quantitative terms, we implement the following equation where �σ max is the maximum stress drop over the rupture area r (which is defined by the strike and dip of the specific finite-fault rupture solution), s(r) is the amount of slip at position r over r , s max is the maximum slip over the rupture area, while û(r) is direction of rake at position r .For each of the finite-fault solutions utilised here 25-28 , we generate 10 5 samples of M e .In doing so, we vary randomly the strike of the surface projection of the rupture by ±3 • (corresponding to a relocation of the surface projection by up to 5 km).We also vary the dip of the fault and the slip rake prescribed by the specific finite-fault solution by ±5 • .We keep �σ max , occurring where slip is maximum, as a free parameter that we vary uniformly in range from 0 to 2 MPa, and from 4 to 6 MPa-in line with global statistics of strike-slip earthquakes 52 .We randomly select values of the controlling parameters within the ranges/sets above, and integrate the infinitesimal torque associated with the infinitesimal stress-drop vector over the rupture area surface.This allows us to generate an ensemble of 10 5 samples of torque variations M e imparted to SC by the Great Wenchuan earthquake rupture.
For calculations of M e that rely on focal-mechanism solutions, we implement the method of Iaffaldano et al. (2022) 15 : the direction of the coseismic stress associated with the occurrence of an earthquake -here indicated with the unit-vector symbol ŝ -can be determined from the strike, dip, and rake (the latter one is a scalar val- ues in the case of focal-mechanism solutions) of the main nodal plane, while its magnitude �σ is kept as a free parameter.The torque variation M e is then parameterised as where r e is the position of the earthquake hypocentre, while e is the size of the rupture area.We constrain the unit vector ŝ from various focal mechanism solutions -i.e., strike, dip, and rake -available in the scientific lit- erature.We elect to assume an uncertainty on each of the parameters equal to 7 • , which is in line with statistical analyses of the focal mechanism parameters of the 2008 Great Wenchuan earthquake 30 .As in the case of calculations relying on finite-fault slip solutions, we take the stress drop �σ variably in range from 0 to 2 MPa, and from 4 to 6 MPa.We draw 10 5 samples of these parameters from uniform distributions within the ranges indicated above, and build an ensemble of 10 5 realisations of M e .We visualise the distribution of M e by plotting the range of sampled magnitudes, assuming two intervals of possible values for �σ max , and the contour containing 99% of the geographical positions where the sampled torque-variation directions �M e /|�M e | intersect Earth's surface -also referred to as torque-variation poles.

W e n c h u a n s e is m ic s e ri e
Figure 1.Velocities observed at continuously-recording (solid arrows) and campaign (empty arrows) GNSS stations used in this study to infer Euler vectors for the South China (SC, in black) microplate motion relative to Eurasia during the periods from January 2002 to December 2004 (in red), and from January 2015 to December 2017 (in blue).All GNSS stations are labeled: italics is used for campaign GNSS stations available during the earlier time period, bold for continuous GNSS stations, and bold-italics for the four continuously recording sites available for both the earlier and later time periods.The black beach ball shows the epicentre and USGS focalmechanism solution of the 2008 M W 7.9 Great Wenchuan earthquake.Turquoise dots are aftershocks.This figure was created utilising the Generic Mapping Tools version 6 (https://www.generic-mapping-tools.org/) 64 .

Figure 2 .
Figure 2. Distributions of velocity change at continuously-recording stations whose data cover both the earlier (from January 2002 to December 2004) and later (from January 2015 to December 2017) periods (a-d), as well as at pairs of closely-located (<30 km) sites (e-h).For the latter ones, data from the first site (campaign) in the label cover the earlier period, while data from the second site (continuously-recording) cover the later period.Distributions are obtained from the difference between ensembles of site velocities sampled from the nominal values and standard deviations of the East and North velocity components.Red dots show the origin (0,0), which would correspond to the centre of the distribution in case of no temporal change of site velocity.Contours show 20%, 68%, and 95% most-recurrent portions of the distributions.See Supplementary Information for maps of sites and site-pairs.

for the period 1 Figure 3 .
Figure 3. (a) Velocity residuals of GNSS stations (empty squares are stations positions) used to infer the SC/EU Euler vector labeled as B3.0P.(b) Same as (a), for the Euler vectors labeled as A3.0.Ellipses are 95% confidence.GNSS stations without residuals have been excluded from the Euler-vector inversion procedure (see main text for details).This figure was created utilising the Generic Mapping Tools version 6 (https://www.genericmapping-tools.org/)64 .
/EU Euler vectors b) Probability that Euler-vector change owes only to data noise

Figure 4 .
Figure 4. (a) Euler vectors for the motion of SC relative to EU before (B3.0P in Table1) and after (A3.0 in Table1) the 2008 Great Wenchuan earthquake, inferred from continuously-recording GNSS velocities (Fig.1).Map shows the positions of the average Euler poles inferred for the two time periods, together with the contours where the most-recurrent 68% (dashed lines) and 95% (solid lines) of the sampled Euler poles fall.Continents are in light grey.SC margins are in dark grey.Inset shows the distributions of sampled angular velocities associated with the Euler poles.(b) Estimate of the probability that random, cross-site correlated noise in GNSS velocities is the only factor determining the difference between inferred Euler vectors for the earlier and later periods.Inset shows the histogram of the obtained F-ratio-test probabilities in noise tests (see main text for details).Part of this figure was created utilising the Generic Mapping Tools version 6 (https://www.genericmapping-tools.org/)64 .

Figure 5 .
Figure 5.Surface velocities calculated at uniformly-spaced, arbitrary positions from the inferred Euler vectors: red arrows are calculated from B3.0P, blue arrows from A3.0 (see Table1).Ellipses around velocity arrows show the 95% confidence.SC microplate margins are in black.Continents are in grey.This figure was created utilising the Generic Mapping Tools version 6 (https:// www.gener ic-mappi ng-tools.org/)64 .

Table 1 .
Euler vectors ( ω ) and associated covariances ( C ij ) for the SC motion relative to Eurasia before and after the 2008 Great Wenchuan earthquake.Labels are as follows: B/A refers to the time period before/after the earthquake.The number following B/A indicates the length in years of the position time-series used to calculate GNSS site velocities.A * appended to the label indicates use of the reduced network (see main text for further details).